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Abstract - Bone remodelling is carried out by 'bone multicellular units' (BMUs) in which active osteoclasts and active 
osteoblasts are spatially and temporally coupled. Bone gain or bone loss from remodelling is directly related to the 
bone balance of individual BMUs. Whilst tetracycline double labelling experiments can reveal part of the dynamics of 
the refilling process taking place in the 'closing cone' of a BMU, less is known about how bone refilling is regulated 
in a BMU. Here, we extend our previous mathematical model of cell development within a single BMU to investigate 
the dynamics of refilling in the BMU and help elucidate the contributions to the refilling rate of osteoblast number 
and osteoblast secretory activity at various positions along the BMU's length. The mathematical model is based on 
biochemical coupling between osteoclasts and osteoblasts of various maturity, and includes the differentiation of 
osteoblasts into osteocytes and bone lining cells, as well as the influence of BMU cavity shrinkage on osteoblast 
development and activity. The model is calibrated against experimental data on cell densities and cell numbers in a 
BMU, and compared with other sources of experimental data for validation, such as tetracycline double labelling data. 
Our model shows that tetracycline double labelling does not reveal the start and the end of the refilling process in 
a BMU. This suggests that at a particular bone site undergoing remodelling, bone formation starts and ends rapidly, 
supporting the hypothesis that osteoblasts behave synchronously. Our model also suggests that part of the observed 
cross-sectional variability in tetracycline data may be due to different bone sites being remodelled by BMUs at different 
stages of their lifetime. The different stages of a BMU's lifetime (such as initiation stage, progression stage, and 
termination stage) depend on whether the cell populations within the BMU are still developing or have reached a 
quasi-steady state while travelling through bone. We find that due to their longer lifespan, active osteoblasts reach 
a quasi-steady distribution more slowly than active osteoclasts. We suggest that this fact may locally enlarge the 
Haversian canal diameter (due to a local lack of osteoblasts compared to osteoclasts) near the BMU's point of origin. 

Key words: bone remodeling, bone multicellular unit, closing cone, matrix apposition rate, tetracycline labeling, 
computational modeling 



1 Introduction 

Bone remodelling renews bone tissue in a spatially and tempo- 
rally discrete fashion by means of 'basic multicellular units' 
(BMUs). In a BMU, osteoclasts (bone-resorbing cells) create 
a resorption cavity called the 'cutting cone' and osteoblasts 
(bone-forming cells) refill this cavity, forming a so-called 
'closing cone' [1,2]. The action of osteoclasts and osteoblasts 
in a BMU is well coordinated such that bone formation closely 
follows bone resorption [1, 3,4]. During remodelling, bone 
may be either gained or lost depending on the balance of 
bone turned over in a BMU, i.e., the amount of bone refilled 
compared to the amount of bone resorbed. Several disorders 
of bone remodelling are associated with bone imbalance in 
BMUs. For example, bone loss in osteoporosis is associated 
with a negative bone balance that is mainly due to an under- 
refilling of the BMU cavity opened up by the osteoclasts [3-6]. 
Bone gain in sclerosing bone disorders such as sclerosteosis, 
van Buchem disease, and high-bone mass phenotype are 
associated with a positive bone balance due to increased bone 
formation [7]. 

The rate of bone formation is determined both by the num- 
ber of active osteoblasts and their level of secretory activity 
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(i.e., the volume of osteoid secreted per unit time by an 
osteoblast) [3, 8]. Whilst bone formation rate in BMUs has 
been extensively studied using histomorphometry techniques 
such as double tetracycline labelling [9-16], much less is 
known of the respective contributions of osteoblast number 
and osteoblast level of activity in a BMU in vivo [17-20]. Un- 
derstanding how these respective contributions are evolving in 
a BMU is necessary to improve our fundamental knowledge of 
the regulation of bone formation, and so of bone balance in a 
BMU. These insights are important for our understanding of 
the development of osteoporosis and for its treatment [5]. 

Bone refilling is known to have a different dynamics in 
cortical BMUs, in which matrix apposition rate decreases 
exponentially in time, than in trabecular BMUs, in which 
matrix apposition rate decreases linearly in time [1, 15,20,21]. 
In Ref. [21], Martin attempted to explain this difference by 
assuming that the different geometries of cortical BMUs (cylin- 
drical) and trabecular BMUs (trench-like or planar) would 
affect the strength of a hypothetic regulatory signal on active 
osteoblasts produced by osteocytes [22]. Martin's model does 
not take into account the number of osteoblasts and their 
level of activity separately. However, it raises the important 
question of the relative importance of different regulatory 
systems: hormones and local regulatory molecules, and geo- 
metric influences on bone formation. The geometry of a BMU 
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may influence bone refilling at several levels and for different 
reasons. It may influence osteoblast secretory activity via 
osteocyte regulation as in Martin's hypothesis [21], and it may 
influence changes in the surface density of osteoblasts due 
to the evolving bone surface area of the refilling BMU cav- 
ity [20,38]. These geometric influences constitute inevitable 
effects that are important to estimate so as to deduce the 
precise contributions of other regulatory mechanisms, such as 
the hormonal and local biochemical regulation of cells. 

The specific influence of geometry on osteoblast devel- 
opment and osteoblast level of activity is hard to determine 
in vivo due to the experimental difficulty of controlling this 
geometry. Geometric changes during new bone formation are 
the result of the (geometry-dependent) collective dynamics 
of osteoblast activity, and also depend on the geometry of 
the bone substrate previously shaped by osteoclasts. In 
Ref. [23], the local curvature of the substrate on which new 
bone is formed was shown to influence matrix apposition 
rate in an in-vitro system, but the respective contributions 
of osteoblast population and osteoblast level of activity was 
not investigated. The difference between matrix apposition 
rate in cortical BMUs and trabecular BMUs also suggests 
that geometry plays an important role during bone apposi- 
tion [1,20,21]. The experimental challenges associated with 
the study of geometric influences on bone cell development 
and on cell secretory activity in vivo emphasise the need for 
accurate computational models to help understand the specific 
contribution of geometry versus other regulatory systems, and 
thereby help interpret experimental data. 

Insights into how osteoblast secretory activity is regulated 
in a BMU could be gained by estimating secretory activity 
levels at different phases of the refilling process in the BMU, 
i.e., at different locations in the closing cone of the pro- 
gressing BMU. The local secretory activity of osteoblasts 
can be deduced from the local measurement of both the 
bone formation rate and the density of cells. Tetracycline 
double labelling and histomorphometric analyses enable the 
determination of geometric properties such as the radius of the 
BMU cavity R and the so-called 'matrix apposition rate' | ^- t R\ 
(i.e., the thickness of the new bone layer deposited per unit 
time, or speed of linear deposition) [1, 16]. The distribution 
of cells along a BMU [24] can be revealed experimentally by 
serial sectioning and histomorphometry [18,25], or confocal 
microscopy [26]. However to our knowledge, this data have 
not been reported in the literature [20]. Marotti et al. [18] 
performed serial sections of BMUs, but a quantitative analysis 
is only presented for three cross-sectional slices. 

In this paper, we therefore follow a different approach 
and take the information of osteoblast distribution along the 
BMU that is predicted using a mathematical model of cell 
development in a cortical BMU. The model is based on 
biochemical regulation of bone cell developments and extends 
our previous model [24] by including important influences on 
the density of osteoblasts in the BMU, namely, a geometrical 
effect of cavity shrinkage on the density of osteoblasts, an 
osteoblast-to-osteocyte transition (that likewise depends on 
the geometry) and an osteoblast-to-bone lining cell transition 
at the end of the refilling process. The model is calibrated 
for osteoclast, osteoblast, and precursor cell numbers in the 
BMU [2, 27], and for osteoblast surface densities at three radii 



of the BMU cavity [18]. The experimental values of osteoblast 
secretory activity kf orm (volume of osteoid secreted per cell per 
unit time) obtained by Marotti et al. [18] at three different radii 
are extrapolated into a functional dependence of kf olm (x,t) 
upon the current radius R(x, t) of the BMU cavity. 

The prediction of our mathematical model for the geometric 
refilling rate are then compared with experimental data from 
tetracycline double labelling experiments [16]. This compari- 
son enables us to analyse data obtained from several BMUs at 
a single snapshot in time [16] in terms of the simulated spatio- 
temporal dynamics of a single BMU. In particular, our model 
suggests that tetracycline double labelling data provides an 
incomplete account of the whole refilling process in a BMU, 
that does not represent the start and the end of refilling. Im- 
portantly, we find that a large variability in experimental data 
may arise if measurements are performed on BMUs at different 
stages of their lifetime, such as initiation stage, progression 
stage, and termination stage (an information that is usually 
not known in histological cross-sections). In particular, the 
different characteristic times associated with the build-up of 
the osteoclast population and the osteoblast population within 
a developing BMU imply that these populations do not reach a 
quasi-steady state at the same time, despite their biochemical 
coupling. This is shown to influence the dynamics of refilling 
at a bone site undergoing remodelling, and so to influence 
tetracycline data collected at this site. 

Below, we begin by providing a general derivation of the 
relationship between cell number, level of secretory activity, 
and geometric refilling rate. This derivation applies to osteons 
of any shape and may thereby be useful for experimental 
studies of the kinetics of regulation of bone balance in BMUs 
(see the Conclusions, Section 5). 

2 Methods 

2.1 Bone refilling dynamics 

Active osteoblasts secrete osteoid, a collagen-rich substance 
that mineralises into new bone. The local production of 
osteoid depends both on the local number of active osteoblasts 
and their level of secretory activity. In the following, we 
consider a transverse cross-sectional slice of thickness 8x, at 
a fixed position x in bone. The slice corresponds to a region 
of interest exhibiting the refilling process taking place in the 
closing cone of the BMU while the BMU progresses through 
bone along the x axis (see Figure la). We denote by S(x,t) 
the surface area of the BMU cavity and by P(x, t) the perimeter 
of the BMU cavity at time t in this cross section. To elucidate 
the relationship between BMU refilling rate, number of active 
osteoblasts, and invidivual osteoid secretion rate, we assume 
that 

1. The osteoid secretion rate kf ovmi (t) (in volume per unit 
time) at time t of a single active osteoblast i in the 
cross-sectional slice only depends on that osteoblast's 
microenvironment (such as the presence of nutrients, 
hormones and regulatory molecules [2, 3, 28], and the 
local curvature of the cavity [20, 23]). This microenvi- 
ronment is assumed to consist of influences from within 
the cross-sectional slice at position x, and so kf 0lmi (t) — 
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Figure 1 - (a) A slice (of thickness 5a:) at position x along the BMU showing the surface area of newly-formed bone, the surface area of the cavity S(x.t) and 
its perimeter P(x.t), and the surface area of the ring-like region containing active osteoblasts S B a (*>?)■ ( D ) Four stages of the formation phase of the BMU: 
1. Onset of formation: differentiation of pre-osteoblasts (OB p ) into active osteoblasts (OB a ); cavity radius equal to the cement line radius R c ; 2. Early stage 
of formation, showing that some osteoblasts have become osteocytes (OCY) embedded in the bone matrix; cavity radius R(x,t); 3. Mid stage of formation; r 
denotes a radial coordinate; 4. End of formation: the active osteoblasts have become bone lining cells (LC) and the cavity radius is the Harversian canal radius 



2. Osteoid is not transported into the cross-sectional slice 
or away from it: an increase in osteoid volume within 
the slice is only due to the active osteoblasts present 
there, and serves entirely to refill the BMU cavity in this 
slice (Figure 1). This hypothesis is supported by the 
observation that osteoblasts do not move significantly 
along the longitudinal axis of the BMU while they deposit 
osteoid [1,2,24,29,30]. 

The total volume of osteoid produced per unit time in the 
slice is given by £,fcf rm i(0> where i = 1)2, ... runs over 
the collection of active osteoblasts in the slice. This newly- 
formed bone progressively closes the BMU cavity, of surface 
area S(x,t), according to: 

j- t S(x,t)8x= - kform, ( (0 = -kformfat) dNo^fat), 
i=l,2,... 

(1) 

where 8No Ba (x,t) is the number of active osteoblasts in the 
slice at x. Introducing the surface density of active osteoblasts, 

8N 0Bi (x,t) 



POB a (M) 



SxP(x,t) 



Eq. (1) becomes: 
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The specialisation of Eq. (3) to a rotation-symmetric cor- 
tical BMU of cavity radius R(x,t) gives, substituting S(x,t) — 
nR(x,t) 2 mdP(x,t) = 2nR(x,t): 
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where \4;R(x,t)\ is the so-called matrix apposition rate [1, 
31]. Interestingly, the BMU cavity radius R(x,t) does not 
appear explicitly in the right hand side of Eq. (4), in contrast 
to the presence of a geometrical property of the cavity (the 
perimeter P) in the right hand side of Eq. (3). In fact, it can 
be shown that the matrix apposition rate vj^ rm at a position r 
of the bone surface, defined by the thickness of newly-formed 
bone deposited per unit time, is always given by 



irrespective of the geometry of the bone surface (such as the 
local curvature) and irrespective of the direction of movement 
of the osteoblasts during the deposition (see Appendix A). 1 
The occurrence of an explicit geometric factor of the BMU cav- 
ity in the right hand side of Eq. (3) reflects the fact that Eq. (3) 
corresponds to the average of the local refilling dynamics 
equation (5) over the cross-sectional slice (see Appendix A). 
In contrast to Eq. (4), Eq. (3) applies to arbitrarily-shaped 
cortical BMUs, as well as trabecular BMUs and formative bone 
modelling for appropriate definitions of S and P. We note here 
that the local law (5) applies in fact more generally to any cell- 
driven biosynthetic or catalytic process occurring from the 
surface of a pre-existing substrate, in which the synthesised 
or catalysed material does not diffuse. 2 

The refilling dynamics (3),(4), and (5) express the general 
relationship between geometric refilling rate (4zS, ^- f R, or 
v form)' osteoblast population 3 (PobJ and osteoblast activity 
(kf orm ). These equations can be used to determine one of these 
quantities from the knowledge of the two others. In particular, 
insights into how osteoblast activity kf orm is regulated in a BMU 
can be gained by determining kf olm at different phases of BMU 
refilling, for example by determining its variation along the 
length of the BMU [18]. Unfortunately, this data has rarely 
been measured [20]. 

Implicit dependences upon geometrical properties of the 
BMU cavity in both osteoid secretion rate kf olm and osteoblast 
density p B ; , may of course be present. In fact, phenomenolog- 
ical relationships between matrix apposition rate and the BMU 
cavity radius in cortical bone [11-15,21], or between matrix 
apposition rate and the local curvature in the in-vitro systems 
of Ref. [23] emphasise such implicit dependences. 



Vf im{r,t) = kotm(r,t)poB a (r,t), 



(5) 



The matrix apposition rate is denoted here by vi as it corresponds to 
the normal component of the velocity of the bone surface, see Appendix A. 

2 For example, Eq. (5) applies to bone resorption if replacing osteoblast 
properties by osteoclast properties. However in this case, the spatial integra- 
tion of the local law (5) does not retrieve an equation of the form (3) due to 
osteoclast movement along the x axis (see also Discussion, Section 4). 

3 Provided that there is no gap between adjacent osteoblasts, the surface 
density of active osteoblasts is the inverse of their so-called 'secretory 
territory', i.e., their contact area with the bone surface, which is measured 
inRefs [17,18]. 
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In this paper, we will investigate more particularly the 
implicit dependence of p OBa upon the BMU cavity radius 
R(x,t) in cortical BMUs by modelling the interactions between 
osteoblasts and the bone surface on the cellular level. The 
implicit dependence of kf orm upon R(x,t) will be extrapolated 
from experimental data from Ref. [18]. 

2.2 Computational model for osteoblast distri- 
bution in the BMU 

The surface density of osteoblasts depends both on biochemi- 
cal processes modifying the number of osteoblasts in the slice, 
and on the local bone surface area, which also evolves in time. 
In the following, we restrict to rotation-symmetric cortical 
BMUs and the refilling dynamics described by Eq. (4). The 
BMU cavity radius is therefore the only geometrical parameter 
to consider. This restriction enables a more direct comparison 
with data derived from tetracycline double labelling experi- 
ments, in which only BMU cavity radii are reported (usually 
averaged over four directions in each cross-sectional BMU, see 
e.g. [9,11]). 

There is an increasing interest in the development of mathe- 
matical and computational models of cell populations in bone 
remodelling, due to the complexity of the spatio-temporal dy- 
namics of cell-cell and cell-bone interactions. In recent years, 
several teams of researchers have developed such mathemat- 
ical models, either focused on the tissue scale [32-37] or on 
the scale of a single BMU [24, 38-43]. To our knowledge, the 
only previous mathematical model investigating the refilling 
process of a BMU in some detail is that of Ref. [38]. However, 
in Ref. [38], the resorption process and the start of the refilling 
process are not modelled. The refilling process starts with an 
assumed initial population of active osteoblasts. 

The present work is based on our previous mathematical 
model of osteoclast and osteoblast development in a cor- 
tical BMU in one spatial dimension (the BMU longitudinal 
axis). In our model [24], biochemical coupling between 
osteoclasts and osteoblasts and migration properties of these 
cells were estimated and led to the emergence of a stable 
moving structure (the BMU). The spatial distributions of cells 
along the longitudinal axis x were segregated, with osteoclasts 
towards the front of the BMU and osteoblasts towards the back 
of the BMU. In the model [24], active osteoblasts derived 
from pre-osteoblasts (OB p ), which themselves derived from 
mesenchymal stem cells, at rates depending on the presence 
of transforming growth factor j3 (TGFp). Active osteoblasts 
were assumed to be eliminated from the 'active' cell pool at 
a constant rate, combining osteoblast-to-osteocyte transition, 
osteoblast-to-lining cell transition and osteoblast apoptosis 
in a single parameter. Here, we extend the model [24] by 
calculating the evolution of the cavity radius according to 
Eq. (4) and by including explicitly the following influences 
on the density of active osteoblasts: 

1. Geometric influence of surface availability: cavity 
shrinkage due to refilling reduces the surface area of the 
bone interface and confines osteoblasts onto a smaller 
surface (or equivalently, into a smaller volume); this 
tends to increase osteoblast density; 

2. Osteoblast-to-osteocyte transition: a fraction of active 



osteoblasts becomes trapped in new bone matrix to 
become osteocytes. This tends to decrease osteoblast 
density; 

3. Osteoblast-to-bone lining cell transition: the end of the 
formation phase at the back of the BMU coincides with the 
terminal differentiation of active osteoblasts into bone 
lining cells. Osteoid secretion is stopped during this 
differentiation; 

4. Osteoblast apoptosis: Due to cavity shrinkage, the vol- 
ume of bone required to deposit a layer of constant 
thickness diminishes as refilling proceeds. At the same 
time, a decreasing number of osteoblasts can access the 
bone surface and remain active. An increasing number 
of osteoblasts that are not active, not becoming osteo- 
cytes nor bone lining cells therefore arises (i.e., 'surplus' 
osteoblasts) [1, 2]. These osteoblasts are believed to 
undergo apoptosis (e.g., due to anoikis). Here a constant 
rate of apoptosis will be assumed for simplicity, however 
with a value calibrated so as to eliminate these surplus 
osteoblasts. 

Following our previous model [24], the material balance 
equation for the volumetric density of active osteoblasts (OB a ) 
in the slice can be written as: 

^OB a (jc,f) =% Bp (x,?)OBp(jc,f)-Ao Bl OB a (x,f) 

+ G(x,t) OB a (x, t ) - <7§?y ■ (x,t). (6) 

In Eq. (6), the source term ^ob p OB p corresponds to the TGFp- 
dependent differentiation rate of pre-osteoblasts (OB p ) into 
active osteoblasts; the sink term — A OBa OB a corresponds to the 
rate of apoptosis (point 4. above). 4 The last two terms in the 
right hand side of Eq. (6) are new compared to the model 
of Ref. [24]. The additional source term GOB a describes 
the tendency of cavity shrinkage to increase the density 
(point 1. above). The additional sink term — (Jq°y describes 
the osteoblast-to-osteocyte transition (point 2. above). The 
precise form of these terms (in particular, their dependence on 
the cavity radius R(x, t )) is presented in detail in Appendix A. 
Details on how the transition of osteoblasts into bone lining 
cells is included in the model (point 3. above) are presented 
below. 

Surface density and volumetric density of active os- 
teoblasts. Biochemical processes of cell development are 
best described in terms of volumetric concentrations of sig- 
nalling molecules and volumetric density of cells. Indeed, bio- 
chemical reaction rates depend on the probability of encounter 
of the interacting biochemical compounds. This probability 
depends in turn on how closely packed the compounds are, 
and so, on their local concentration or density (number per 
unit volume) [44]. 5 

4 Whilst the sink term — A OBll OBj has the same form as in Ref. [24], it 
now solely represents apoptosis of active osteoblasts, and not other depletion 
pathways of the pool of active osteoblasts (such as the osteoblast-to-osteocyte 
and osteoblast-to-bone lining cell transitions). 

5 To align with common practice, we use the terminology 'concentration' 
for signalling molecules and 'density' for cells, even though both terminolo- 
gies refer to the same units (i.e., number per unit volume). 
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Table 1 - Osteoid secretion rate vs BMU cavity radius (in dogs) 



R d °£ [um] 


*£L [MmVday] 


15 





17 


90 


31 


103 


65 


180 



The material-balance equation (6) for active osteoblasts 
involves biochemical reactions depending on volumetric 
cell densities (e.g., OB p ) and molecule concentrations (e.g., 
TGFp) [24]. This equation is therefore written for the vol- 
umetric density of active osteoblasts, OB a . However, active 
osteoblasts in a BMU are only found on the bone surface and 
so the refilling dynamics described by Eqs (3)-(5) involves 
their surface density Pob.,- To relate the volumetric density 
OB a to the surface density p OBa such that OB a conveys the 
information of how densely packed active osteoblasts are on 
the bone surface, we use the fact that active osteoblasts in a 
BMU form a single layer of cells against the cavity walls [18]. 
In the cross-section, active osteoblasts are confined to an 
annulus S OB l of width /i OBii corresponding to the typical height 
of an active osteoblast (see Figure la). Within this annulus 
SoB a 7 active osteoblasts are distributed fairly uniformly, and so 
SoB a constitutes an appropriate representative surface element 
to define the volumetric density of active osteoblasts [44, 
Chap. XIV], i.e.: 6 



8N 0Bi (x,t) , P(x,t) 

OB a (x,t) = '- = p 0Ba (-Mh 



8x S, 



OB a 



(7) 



where Eq. (2) has been used for the second equality. In a 
rotation-symmetric BMU, P(x,t) = 2%R(x,i) and 



S OBi Xx,t) = n[R(x,t) 2 - {R(x,t)-h OB J 2 ] . 



(8) 



The mathematical relationship between osteoblast surface 
density and volumetric density can thus be written as: 



PoBaCV) =OB a (*,f) ^OB a 1 - 



<OB a 



2R(x,t) 



(9) 



We note that in the above relationship between p B a and OB a , 
the geometric factor (l — 2R(xt) ) accounts f° r the nonzero 
curvature l/R of the cavity (a flat surface has zero curvature, 
corresponding to the limiting case R °°). Substituting 
Eq. (9) into Eq. (4) gives: 



f t R{x,t) = -k form (x,t) OB a (x,f) /z 0Ba 1 



'OB a 



2R(x,t) 



(10) 



Osteoblast-to-bone lining cell transition and osteoid se- 
cretion rate. At the end of the formation phase at the 
back of a BMU, active osteoblasts terminally differentiate 
into bone lining cells, thus forming a cellular layer on the 

6 Choosing for example the cavity surface area S(x,t) or the constant 
osteon cross-section area rather than S B a to define the local volumetric 
density OB a (x,r) in Eq. (7), would not make OB a (.v,r) represent how densely 
osteoblasts are packed on the bone surface, an information that is impor- 
tantly conveyed by the surface density p B a in Eqs (3)-(5). Indeed, active 
osteoblasts are not distributed uniformly across these other representative 
surfaces. 
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10 20 



30 40 50 60 70 
cavity radius R [um] 



80 90 100 



Figure 2 - Assumed functional relationship kf 0Tm (J?) bewteen osteoid secre- 
tion rate k[ orm and BMU cavity radius R (solid line). Data points (crosses) are 
taken from Table 1 and scaled according to Eqs (11), (12). 

newly formed bone surface [1, 31]. To account for such a 
transition in the model while ensuring the spatial exclusion 
of active osteoblasts and bone lining cells, we assume that 
active osteoblasts whose osteoid secretion rate kf 0lm is rapidly 
dropping become bone lining cells. Osteoid secretion rate by 
individual osteoblasts varies as refilling proceeds in cortical 
BMUs [18, 19]. The precise factors that influence osteoid 
secretion rate are currently not well known. The secretion rate 
depends on the amount of organelles involved in glycoprotein 
and protein synthesis, and so is likely to depend on the cell's 
protoplasmic volume [18, 19,25]. The observed flattening of 
osteoblasts as refilling proceeds [1, 18, 19, 31] is associated 
with a reduction in cell volume and in osteoid secretion 
rate [18,19,25]. 

To account for the variation of osteoid secretion rate during 
the refilling process in our model, we assume that kf orm 
depends on the current radius of the cavity: kf orm (x,t) — 
kf 0lm [R(x,t)) , with kf olm (Rn) = 0, where R^ is the Haver- 
sian canal radius at which osteoblasts have become (non- 
synthesising) bone lining cells (see Figure lb). Marotti et 
al. [18] have reported values for the osteoid secretion rate at 
three different radii in canine cortical BMUs. We add to this 
data a zero secretion rate at a Harversian canal radius /?H dog ~ 
15 um, corresponding to the average Haversian canal radius 
in dogs [31] (see Table 1). The canine data of Table 1 is then 
used to generate pseudo-human data by appropriate scalings 
as explained in the next paragraph. We finally interpolate 
this scaled data by a continuous piecewise polynomial curve 
kf olm (R) that has a sharp stepwise increase from R — R^ = 
20 um to R = 22 um, and a linear slope for R > 22 um, as 
depicted in Figure 2. 

Scaling data from animal models to humans. 

Tetracycline-based bone histomorphometry has been 
performed in several species, such as humans [45, 46], 
sheep [16], dogs [11, 13-15, 18], and cats [12]. However, 
some of the quantitative measurements of BMU-related 
quantities useful to define, calibrate and validate our model 
are derived from animal models, such as cell surface densities 
and osteoid secretion rates at different BMU cavity radii [18], 
total cell numbers [27,31], and tetracycline double labelling 
radii [16]. To utilize such data consistently with human 
data in the calibration of the model to a human cortical 
BMU, we scale some of this animal data to (pseudo-)human 
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values according to known (or suspected) cross-species 
differences [47]. 

In particular, we scale the radius data R d ° s from Table 1 to 
(pseudo-)human radius data R by assuming a linear relation- 
ship # = R(R d0 %) such that R(R H dog ) =Rh mdR(R c do &) =R C , 
where fl H dog » 15 um, Rft rj 20 um are typical Haversian 
canal radii in dogs and humans, and /? c dog « 70 um, R c w 
100 um are typical cement line radii in dogs and humans, 
respectively [1,31,47] (see Figure lb). 7 This gives: 



R(R Ao Z) =Rn+ ,f c (R d °z - V og ) ■ 



(11) 



We also scale the osteoid secretion rates k d °^ of Table 1 
estimated by Marotti et al. [18], by a factor ttk !orm = 1.25 to 
account for higher secretion rates in humans [38], i.e.: 



fcform — ^foim^form- 



(12) 



The specific numerical value of the factor Ct^ lbmi is chosen so as 
to reduce the total number of active osteoblasts in the BMU (see 
Results and Discussion, Sections 3,4). The scaled radius and 
osteoid secretion rate data of Table 1 are shown in Figure 2 
(crosses). 8 

Scalings are also applied to the data of Table 2 reporting 
canine osteoblast surface densities p^ at three BMU cavity 
radii (see Results, Section 3). The radii R dog are scaled 
to correspond to human average BMUs as in Eq. (11). The 
osteoblast surface densities p^ 8 are scaled to pseudo-human 
data as p B a = o;^ 1 Po^f- The choice of the scaling factor 
ar 1 = 0.8 is motivated by the corresponding scaling factor 

Morm 

°Vm = 1 - 25 performed on kf ^ m . Indeed, increased osteoid 
secretion rate is associated with increased cell volume, and 
thus with a corresponding decrease in cell density [18, 19]. 

Governing equations. The equations describing the evolu- 
tion of the BMU cavity radius (10) and the evolution of the 
density of active osteoblasts (6) both govern the dynamics 
of the refilling towards the back of the BMU. Following 
our model [24], we define three other state variables of the 
system: the density of pre-osteoblasts OB p (x,f), the density 
of active osteoclasts nuclei OC a (x,f) and the concentration 
of transforming growth factor |3, TGFp(x,f). 9 The balance 
equations governing these quantities involve several addi- 
tional signalling molecules, i.e. the receptor-activator nuclear 
factor >cB system (RANK, RANKL, and OPG) and parathyroid 
hormone (PTH). However, the dynamics of these additional 
signalling molecules is fast compared to the characteristic 
times of cell behaviours [24]. Their balance equations can 
be taken in a quasi-steady state, enabling the concentration of 
these signalling molecules to become algebraic expressions of 
the variables OB p ,OB a ,OC a and TGFp. The material balance 
equations governing the evolution of all the variables are 
presented in Appendix A, see Eqs (19)-(30). 

7 In Ref. [31], the average cement line radius in dogs is reported to be 
« c d °s S3 60 um. In Ref. [47], it is reported to be K c do « ss 77 um. Since one of 
the three BMU cavities in Ref. [18] has a radius 65 um in the formation phase, 
we take here 70 um as the representative cement line radius. 

8 The model of Polig and Jee [38] suggests that may take values 

from 1.27 to 1.7. 

'Here and in Ref. [24], OC a s represent mononucleated entities incorpo- 
rated in a multinucleated active osteoclast. Multinucleated active osteoclasts 
are composed of about 10 nuclei [27]. 



Table 2 - Osteoblast surface density vs BMU cavity radius (in dogs) 

fl d °g [um] pp°f [mm' 2 ] 

~. 17 4500 

31 8770 

65 10000 



3 Results 

Numerical simulations. The material balance equations of 
the average densities of cells and concentrations of signalling 
molecules (19)-(30) (see Appendix A) together with Eq. (10) 
form a coupled system of five partial differential equations 
(PDEs). These PDEs were solved numerically to evolve the 
system for 150 days from an initial condition with a small, 
localised population of active osteoclast nuclei (OC a s) and 
given distributions of OB u s and OC p s concentrated around the 
(growing) blood vessel extremity of the BMU, as in Ref. [24], 
see Figure 3a. The initial cavity radius was set to the cement 
line radius R c as only the refilling of the cavity is considered 
in this paper. We assumed a uniform density of osteocytes 
OCY exp .(x,r) = 20000/mm 3 [31]. 

The cells present initially at t = generate a biochemical 
positive feedback loop, captured in the governing equations: 
the OC a s free TGFp from the bone matrix into the microen- 
vironment, which promotes the differentiation of OB u s into 
OB p s. The ligand RANKL expressed on OB p s promotes in 
turn the differentiation of OC p s into OC a s, thus enabling the 
process to be sustained. Further towards the back of the BMU, 
the concentration of TGFp drops [24]. This facilitates the 
differentiation of OB p s into OB a s [2, 3]. At this point, bone 
formation starts and the cavity refills. After an initial transient 
during which these processes take place and the population of 
osteoblasts increases (see Figure 3b), the cell densities reach 
steady-state distribution profiles progressing forward in bone 
without changing shape (Figure 3c). Meanwhile, the bone 
cavity progressively refills (as shown at the bottom of each 
graphs in Figure 3a-c). The refilling rate is strongly reduced 
when the cavity radius reaches values close to the 'target' 
Haversian canal radius R^ = 20 um. Indeed, close to R^, 
the osteoid secretion rate kf olm drops according to Figure 2. 
The distribution of kf orm (R(x,t)) along the BMU's longitudinal 
axis x at t = 150 days is shown in Figure 3d, along with 
the locations at which k{ olm is 90% and 10% of the step 
height in Figure 2. These locations suggest the transition of 
active osteoblasts (OB a s) into bone lining cells (LC) and are 
represented in Figure 3c by the start and end of the OB a — » LC 
transition arrow. The dot-dashed line in Figure 3c repre- 
sents the active osteoblast surface density p B a and constant 
bone lining cell density p LC = 2300/mm 2 (see right vertical 
axis) [31]. Surface densities p OBa < 2300/mm 2 are shaded 
as they do not correspond to real active osteoblasts: at these 
locations, osteoblasts are assumed to have become bone lining 
cells of constant density no longer secreting osteoid and no 
longer undergoing apoptosis. 

Calibration of the cell distribution profiles. Quantitative 
data on the distribution of cells in cortical BMUs, and/or 
their total number, is relatively sparse [2, 18, 20, 27, 48]. We 
calibrated our model such that the cell distribution profiles are 
in reasonable agreement with the following: 
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Figure 3 - (a)-(c) The profiles of cell densities along the BMU's length (top of 
each panel) and the corresponding shape of the closing cone (bottom of each 
panel) are shown at time (a) t = days (initial condition), (b) ? = 75 days, and 
(c) f = 150 days. To ensure that all cell densities are visible, OBu, OC p , and 
OC a were multiplied by a factor 20. The active osteoblast surface density p B u 
and bone lining cell surface density p LC are also shown (dot-dashed lines). In 

(b) , individual osteoid secretion rate kf mm has not yet dropped enough for 
OB a to differentiate into LC. In (c), an OB a — > LC transition takes place as 
indicated by the arrow. The arrow starts and ends where kf <xm (x,t) is 90% 
and 10% of the "step height" in Figure 2 (see also panel (d)). The crosses in 

(c) correspond to the experimental data of osteoblast surface density from [18] 
listed in Table 2, scaled to human values, (d) Profile of osteoid secretion rate 
kf orm (x,t) by active osteoblasts at t = 150 days. The grey lines indicate where 
harm is 90% and 10% of the "step height" in Figure 2 (OB a — > LC transition). 



• The total number of precursor cells (N b u , N 0Cp ) in steady 
state is about 20; 

• The total number of active osteoclast nuclei (N 0C J in 
steady state is about 100 [27]; 

• The surface density of active osteoblasts (PoB a ) in steady 
state coincides with data reported by Marotti et al. [18] 
at three different cavity radii, whilst the total number of 
active osteoblasts (N OB J in steady state is in the range 
2000-6000 [27,31,38] (see also Discussion, Section 4). 

The osteoblast surface densities p^ 8 reported by Marotti et 
al. [18] in three canine BMU cross-sections are shown in 
Table 2. The radii R Ao ° and surface densities p B a were scaled 
to correspond to human average BMUs as mentioned above. 

The data on osteoblast surface density at different cavity 
radii thus rescaled is shown in Figure 3c (crosses) along with 
the calculated surface density profile p OBa (blue dot-dashed 
line). 

To calibrate the model against the number of active osteo- 
clast nuclei and precursor cells in the steady-state BMU, the 
cell densities OC a (x,f), OC p (x,f ) and OB u (jc,f ) were integrated 
over the BMU cavity space R(x,t) at t = 150 days. An 
appropriate rescaling of model parameters was thereby deter- 
mined (see Appendix A) and the model parameters modified 
according to Eqs (32). With the model parameters listed in 
Table 3 in Appendix A, we obtained: 



oc a 



100, No 



: 20, No 



20 



(13) 



To calibrate the model against the active osteoblast surface 
densities at three different cavity radii (Table 2) and against 
the total number of active osteoblasts, we modified the os- 
teoblast apoptosis rate A B a , the model parameters related to 
the scaling factor Ct B a , an d tne (canine-to-human) scaling 
factor 0f£ form . The apoptosis rate A B a influences the spatial 
rate of decrease of the osteoblast distribution at the back of 
the BMU [24]. The factor a B a scales the value of the density 
at each location along the BMU uniformly. The scaling factor 
(Xk folm enables to modify the total number of active osteoblasts 
in the steady state BMU (e.g., increasing OLk {oTm reduces N OBa ). 
The data could be matched very well by modifying these two 
properties of the distribution profile, as seen in Figure 3c. 
The total number of pre-osteoblasts (obtained by integrating 
OB p (jc,f) over the cavity space) and the total number of 
active osteoblasts (obtained by integrating the surface density 
p OBli (x, t) over the cavity-bone surface) in the steady state 
BMU (at t — 150 days) were: 



No 



11100, N 



■ 5200. 



(14) 



The calibrated model was then validated against three in- 
dependent sets of experimental data: (i) the density of bone 
lining cells reached at the end of the formation phase at the 
back of the BMU; (ii) the length of the formation cone of the 
BMU; (iii) tetracycling double labelling data. 

Osteoblast-to-bone lining cell transition. In the model, 
refilling is strongly reduced when the radius of the cavity 
reaches the Haversian canal radius (Figure 2). The surface 
density of active osteoblasts at which this occurs should be 



7 



— 

Pi 

< 



1.5 



0.5 




Exp. data [16] (scaled) 
x = -3.5 mm (exact MAR) 
x = —2.5 mm 
x = —3.5 mm 
x = —4.5 mm 
x = —4.65 mm . * 
x = —4.7 mm / 



20 40 60 

cavity radius R\ [urn] 



80 



100 



Figure 4 - Matrix apposition rate versus BMU cavity radius obtained from 
our model (lines) and from tetracycline experiments (dots). The exact matrix 
apposition rate | §jR(x, t) | versus R\ = R(x, t) is shown at position x = —4 mm 
in bone (solid grey line). The approximate matrix apposition rates | Rl ^ t R[ | = 
«(.r,t+At) R{x.t) | yersus ^ _ R^x^t) are shown at .v = —2.5 mm, —4 mm, 
—4.5 mm, —4.62, and —4.7 mm for a fixed interval A( = 10 days. Data 
points are taken from Ref. [16, Fig. 2] and have been appropriately rescaled 
to human average BMU sizes (see text). 

consistent with experimentally observed surface density of 
bone lining cells, p LC ~ 2300/mm 2 [31]. We see from 
Figure 3c that this constraint was satisfied by the model. 
Indeed, the surface density p 0Ba reached the value 2300/mm 2 
precisely in the zone in which osteoid secretion rate dropped 
from 90% to 10% of the step height of k form (R) around RxiRu 
in Figure 2 (see also Figure 3d). 

Length of the formation cone. Considering the start of the 
formation cone in Figure 3c to be around x w —0.5 mm and 
the end of the formation cone to correspond to the intersection 
of the active osteoblast surface density and bone lining cell 
density at around x w —3 mm, the length of the formation 
cone is about 2.5 mm. This corresponds to values estimated 
for humans [31]. 

Tetracycline double labelling data. Tetracycline double 
labelling is an experimental technique enabling the estimation 
of the matrix apposition rate ^R. Tetracycline incorporates 
actively mineralising tissue with a fluorescent dye [1,9,31]. 
Its injection at two successive time points t\ and ?2 reveals the 
radii of the BMU cavity at these time points, R\ and R2, from 
histological slices extracted later [11, 12, 16]. The average 
matrix apposition rate between t\ and ?2 can be estimated as 
MAR = J R t 2 2 Zfi I ■ The data is usually presented in either of 
two forms: a plot of the data pairs (R\,R2) [12, 15, 16] or a 
plot of MAR vs either R 1 or R 2 [11, 13, 14, 16]. In both cases 
a linear relationship is experimentally found. 

In Figure 4 we plot the 'exact', instantaneous matrix appo- 
sition rate | §~ t R(x, t ) | versus R(x, t) at x = —4 mm (solid grey 
line) and the approximate, average matrix apposition rates 
I R ( x ' t+At 2~ R< " x '^ I versus R(x,t) obtained in our simulations 
during the BMU cavity refilling when t increases from to 
150 days. The effect of sampling the radii at various positions 
x in bone is shown (black lines), whilst the time interval 
At = 10 days (corresponding to the time interval between 
the tetracycline injections) is kept fixed. These results are 
compared to experimental data from Metz et al. [16]. This 
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Figure 5 - Parametric plots of cavity radius pairs (i?i,i?2) with R\ = R(x,t) 
and i?2 = R(x, t + Af) obtained from our model at a fixed location in bone 
(x = —4 mm) when t evolves from to 150 days. The figure shows the 
effect of choosing different time intervals A; = 10 days, 12 days, and 15 days 
(corresponding to typical time intervals between tetracyline injections). The 
calculated curves intersect with the diagonal at the onset of formation and at 
completion of formation. 

data was collected in sheep and plotted in Ref. [16, Fig. 2] as 
MAR versus "percent unfilled", i.e., the percent of bone thick- 
ness remaining to deposit so as to reach the sheep Haversian 
canal radius. In Figure 4, we rescaled this "percent unfilled" 
data to human-sized BMU cavity radii similarly to Eq. (11), 
such that 100% unfilled corresponds to the cement line radius 
R c = 100 um and 0% unfilled corresponds to the Haversian 
canal radius Rn = 20 |jm. 

In Figure 5 we plot the trajectory of the radii pairs 
(Ri,R2) = (R{x,t),R(x,t + At)) obtained during the BMU 
cavity refilling at x = —4 mm when t increases from to 
150 days for various time intervals At. 



4 Discussion 

Cell distribution profiles within the BMU. As in our pre- 
vious model [24], the distribution of cell densities along the 
longitudinal spatial coordinate x evolves from the initial con- 
dition into a stable multicellular travelling-wave-like structure 
(the BMU) progressing through bone without changing shape, 
with osteoclasts towards the front of the BMU and osteoblasts 
towards the back of the BMU (Figure 3a-c). This developing 
and progressing structure models the initial and quasi-steady 
phases of a cortical BMU. Its emergence is due to the biochem- 
ical coupling between the bone cells and the cells migration 
properties [24]. 

The most important difference in the cell distribution pro- 
files of active osteoclasts (OC a ) and pre-osteoblasts (OB p ) in 
Figure 3a-c compared with Ref. [24] is a vertical scaling 
factor due to the calibration of total cell numbers performed 
here. The overall shape and spatial extension of the distribu- 
tions of OC a s and OB p s is otherwise unaltered. By contrast, 
the shape of the cell distribution profile of active osteoblasts 
(OB a ) is appreciably modified compared to Ref. [24, Fig. 2] 
due to the addition of the geometric influence of surface 
availability, osteoblast-to-osteocyte transition, and osteoblast- 
to-bone lining cell transition in the model. In the current 
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model, active osteoblasts rise more sharply to their maximum 
value at the start of formation (near x = —0.5 mm). The 
decrease of the density of OB a s towards the back of the BMU is 
also more linear in Figure 3c compared to the exponentional- 
like decrease in Ref. [24, Fig. 2], an effect attributable to 
the tendency of cavity shrinkage to concentrate the density 
of active osteoblasts. The generation of osteocytes also has 
an important influence on the decrease of the density of OB a s 
towards the back of the BMU. Indeed, one can observe a small 
transient reduction in the decrease of OB a as a slight "bulge" 
in the OB a density profile near the OB a — > LC transition, where 
osteocytes stop being generated. This "bulge" is generated 
despite the density-concentrating effect of cavity shrinkage 
being reduced at the same time. 

The calibration of the active osteoblast surface density 
PoB a (dot-dashed blue line) to Marotti et a/.'s data (crosses) 
in Figure 3c could be performed very well, in particular 
thanks to the pseudo-linear decrease of the active osteoblast 
population towards the back of the BMU obtained with the 
model. Remarkably, this calibration is consistent with the fact 
that the cavity is refilled up to the Haversian canal radius Ra 
precisely when the surface density p OBa reaches a bone lining 
cell surface density of 2300/mm 2 [31]. This is indicated in 
Figure 3c by the fact that the p OBa curve intersects the OB a — > 
LC transition arrow. With poorly-calibrated model parameters 
the cavity was either under-refilled when p 0Ba reached the 
value 2300/mm 2 or refilled too early, at values of p OBii larger 
than 2300/mm 2 , leading to lengths of the closing cone shorter 
than reported values (see also below). 

Osteoblast apoptosis in a BMU is believed to serve to elim- 
inate 'surplus' osteoblasts, i.e., osteoblasts that are neither 
forming bone, nor becoming osteocytes, nor becoming bone 
lining cells [2]. The rate of osteoblast apoptosis, as a reg- 
ulatory mechanism of osteoblast number in the BMU closing 
cone, is thus likely to depend on the phase of the refilling 
process in the BMU. Indeed, as the cavity shrinks, the number 
of osteocytes to be generated per unit time decreases and 
osteoblasts are confined in a smaller volume, resulting in an 
increase in surplus osteoblasts. Nevertheless, the constant rate 
of osteoblast apoptosis A 0Ba assumed in our model enabled 
us to obtain a profile of active osteoblast surface density 
that accurately matched the experimental measurements of 
Marotti et al. [18] and the bone lining cell density. 

Number of osteoblasts within the BMU Precise estimates 
of the total number of active osteoblasts in a BMU are difficult 
to find in the literature. Numbers in the range 2000-4000 [27], 
of about 3000 per millimetre of BMU length [38], and up to 
6750 [31, Table 6] are reported, probably reflecting cross- 
sectional and/or cross-species variabilities. Estimates of the 
local surface density of active osteoblasts PoB a , or of their 
secretory territory (corresponding to pQ B if closely-packed) 
are likely to be more reliable as they are based on local 
measurements only and do not depend on the length of the 
closing cone of the BMU. A lower bound of the total number 
of active osteoblasts A^,B a can be estimated based on a lower- 
bound surface density of 2300/mm 2 (surface density of bone 
lining cells [31]), and a lower-bound cavity-bone surface area 
given by a linearly changing cavity radius along the formation 
cone of the BMU (straight lines being the shortest). With a 



formation cone length Lf orm ~ 2.5 mm [31], Haversian canal 
radius Rn ss 20 um and cement line radius R c w 100 um, this 
provides the lower bound Af B a ^ 2170. Following Parfitt [31], 
a reasonable upper-bound estimate may be found similarly 
by assuming a maximum cavity-bone area given by that of a 
cylinder of radius R c and length Lf orm and an average surface 
density p OB ., ~ 4500/mm 2 . This provides the upper bound 
Nob.„ < 7000. 

The total number of active osteoblasts in the steady-state 
BMU obtained in our simulations (N OBa ~ 5200) falls well in 
the above ranges. This total number of active osteoblasts 
was reached for a scaling of the canine experimental osteoid 
secretion rates and osteoblast surface densities by the factors 
0Ck fmm — 1-25 and a^ 1 — 0.8, respectively, as explained in 
Methods, Section 2. Without these scalings (for a BMU 
properly calibrated to the other data), the number of active 
osteoblasts was about 6600. 

As mentioned above, the length of the BMU closing cone 
is an important factor determining the total number of active 
osteoblasts. The BMU progression rate in bone, v BMU , was 
found to directly affect this length, and so also the total 
number of active osteoblasts. Values of v BMU in the range 20- 
40 um/day [2, 31] are cited in the literature. The simulations 
presented here were performed with a BMU progression rate 
Vbmu = 30 um/day. With a BMU progression rate v BM u = 
40 um/day and no scaling of the canine data {i.e., Ofjfc fomi = 1), 
about 8800 active osteoblasts were found in steady-state in 
a closing cone measuring about 3.5 mm. These results show 
that a great variability in the total number of active osteoblasts 
can be achieved within the model with still physiologically 
reasonable parameter values. Such variability can therefore 
also be expected in experimental measurements. 

Finally, the total number of pre-osteoblasts in the BMU in 
our simulations (N OBp ~ 1 1 100) may seem high, although we 
are not aware of pre-osteoblast cell counts in BMUs reported in 
the literature. The high value obtained for N OBp compared to 
N OBsi is mainly due to the integration of the volumetric density 
OB p across the whole cavity surface S(x,t), as opposed to the 
integration of the density of active osteoblasts near the cavity- 
bone surface only. We note, however, that the fate of all 
pre-osteoblasts in the model is to become active osteoblasts 
at some point of the refilling process, i.e. pre-osteoblasts are 
not assumed to undergo apoptosis. In fact, the number of pre- 
osteoblasts in a cross-sectional slice of thickness 8x = h OBlL = 
15 um at a fixed position in bone reaches a maximum value of 
about 390 near the start of the refilling process, then quickly 
decreases to zero as pre-osteoblasts become active. This 
activation process forms a layer of active osteoblasts against 
the bone surface with precisely the correct surface density, 
as indicated by the match of the surface density of active 
osteoblasts p B a with the experimental data in Figure 3c. 
This shows that in terms of cell densities, the model is well 
calibrated, but that absolute number of cells derived from the 
model by integration of the densities may be enticed with large 
uncertainties, associated with uncertainties (or variabilities) in 
the domains of integration. 

Shape of the BMU cavity (closing cone). The bmu cavity 
represented below each plot in Figure 3a-c has to be in- 
terpreted with some care as only the refilling process was 
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accounted for, not the resorption process. A major difference 
between the resorption process by active osteoclasts and the 
formation process by active osteoblasts within a cortical BMU, 
is that osteoclasts move longitudinally (along x) while they 
resorb bone [31]. This complicates the governing equation for 
the cavity radius in the cutting cone. Indeed, transport terms 
along x due to osteoclast movement are present, which involve 
spatial derivatives. We note that this complication does 
not necessarily arise in so-called 'type II osteons' (formed 
by BMUs following a pre-existing Haversian canal [31,49]) 
and in trabecular remodelling as in these situations, cellular 
resorption may be progressing radially only [31]. 10 

The local shape of the Haversian canal resulting from the 
passage of a BMU in a region of bone undergoing remodelling 
depends on the stage of the BMU's life that this region of 
bone experiences. The extent of both resorption and formation 
depends on whether the BMU was in an initiation phase (early 
life), in a quasi-steady state (mid life) or in a termination phase 
(late life). In fact, several levels of quasi-steady states can be 
defined in a BMU. Indeed, whilst the population of certain cell 
types may have reached a quasi-steady-state within the BMU, 
the population of other cell types may still be developing. Due 
to the longer lifespan of osteoblasts compared to osteoclasts, 
we observe that osteoclast cell densities reach a steady state 
earlier in the BMU's life than osteoblast cell densities. For 
example, in Figure 3b (at t = 75 days), the density of active 
osteoclasts (OC a ) has already reached a quasi-steady state but 
not the density of OB a s near the back of the BMU. It is therefore 
likely that regions of bone in a BMU's path experience the full 
extent of the BMU's resorption but are only partially refilled 
due to the density of osteoblasts having not fully evolved to 
a steady state at this location. This situation is mimicked 
in our simulations by the refilling process occurring near 
x = —4.65 mm in Figure 3a-c, where a large cavity generated 
by a mature BMU osteoclastic front is only partially refilled. 1 1 
This creates a local enlargement of the Haversian canal that 
may be contributing to local spatial nonuniformities observed 
in Haversian canal sizes [43,50,51]. The observation of such 
local enlargements may be useful for the determination of the 
location of a BMU's birth from micro-CT scans exhibiting the 
Haversian pore system [51]. 

Comparison with tetracycline double labelling data. The 

calibration of the active osteoblast surface density p OBa (dot- 
dashed blue line) to Marotti et a/.'s data (crosses) in Fig- 
ure 3c also provides a remarkable consistency between the 
BMU cavity refilling dynamics of the model and experimental 
tetracycline double labelling data (see Figure 4). However, 
the prediction of the model and the experimental data need 
to be compared with some care. The instantaneous matrix 
apposition rate \4zR(x,t)\ (Figure 4, solid grey line) differs 



For example, in the mathematical models of trabecular BMUs of Refs [39, 
40,42], only radial resorption is considered. 

1 1 In the simulation corresponding to Figure 3a-c, the density of osteoclasts 
quickly increases from almost zero to stabilise to a quasi-steady state at 
around 7-8 days (data not shown). All the bone in the BMU's path located 
at x > —4.65 mm therefore experiences the full extent of the resorption phase 
of the BMU, where —4.65 mm corresponds to the location of BMU origination 
xq = —4.85 mm plus the distance travelled by the BMU during 7 days (about 
200 pm). In contrast, the density of osteoblasts takes about 50 days from the 
BMU's origination to increase to its maximum value. 



significantly from its approximation | ^w+^j-^w) | (solid 
black lines and dots). Experimentally, relatively large time 
intervals between the tetracycline injections (Af w 7-12 days) 
are usually prescribed so as to clearly distinguish both labels 
in the histological sections [10]. However, mathematically 
| J^/?(x,f)| is defined in a continuous limit in which At — > 0. In 
fact, the experimental matrix apposition rate above represents 
a better estimate of \§jR(x,t + y)| rather than \^R(x,t)\, 
but it is nevertheless usually compared with the radius of the 
cavity at time t . 

Experimental measurements of label radii are often per- 
formed on a collection of different osteons, rather than on 
serial sections of a single osteon [1 1, 12, 16]. Having observed 
great cross-sectional variations of cement line radii and Haver- 
sian canal radii, the authors of Ref. [16] have presented MAR 
versus the percentage of unfilled bone to normalise the data 
from different osteons. We therefore took this normalised data 
to compare with the result of our simulations performed on a 
single BMU. 

The numerical results in Figure 4 (lines) clearly exhibit the 
start and the end of the formation phase, at which MAR = 0. 
The sharp rise of the active osteoblast population at the start 
of formation is reflected by a sharp increase of MAR near 
the cement line radius R c = 100 |jm. Matrix apposition rate 
then steadily decreases as refilling proceeds and cavity radius 
shrinks. The Haversian canal radius eventually reached is 
given by the point at which MAR = 0. Because osteoid 
secretion rate is assumed to decrease sharply when the cavity 
radius approaches R^ = 20 |jm, a sharp fall of MAR near 
Rii is also observed. The rise of MAR at the start of the 
formation and the fall of MAR at the end of the formation are 
not reflected in the experimental data. This strongly suggests 
that both the start and the end of the formation phase may 
occur rapidly, perhaps as a synchronised event amonsts all 
the osteoblasts. A synchronised halt of formation may allow 
the BMU cavity to reach a target size more consistently. One 
potential mechanism is that osteoblasts sense the proximity of 
the blood vessel or other localised components, as suggested 
e.g. in Ref. [20]. 

Interestingly, for the most part, the numerical curve MAR 
versus Ri exhibits a pseudo-linear character. This pseudo- 
linear character is consistent with the experimental literature, 
in which such apparent linear relationship is often fitted by 
lines passing at or near the origin (R, MAR) = (0,0) or 
is written as the phenomenological law: MAR = jnR = 
CR [11,12, 14, 15,38], despite MAR being normally zero at the 
Haversian canal radius Rh (not at R = 0). This gives additional 
support to the suggestion that active osteoblasts may stop bone 
formation in unison. Indeed, the continuation of the pseudo- 
linear part of the solid black curve in Figure 4 reaches a point 
near the origin, while the end of formation clearly corresponds 
to the point (R, MAR) = (7?h,0). This thus suggests that near 
the end of formation, the curve MAR(7?) falls abruptly to 
the point (/?h,0), as in Figure 4. Similar remarks apply to 
the (R\,R2) data in Figure 5. In Figure 5, the start and end 
of formation correspond to the points (R C ,R C ) and (Rh,Rh), 
respectively. A pseudo-linear behaviour is also observed for a 
significant part of the numerical curves, whose extrapolation 
may pass at or near the origin (R1.R2) = (0,0). 

In spite of the normalisation of the data performed by 
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Metz et al. [16] to reduce cross-sectional variability and 
thereby approximate a 'typical' osteon, another source of 
variability in the data is manifestly present. Our numerical 
model suggests that the stage of the BMU's life at which the 
data is collected may be an important contributor to such 
variability. In Figure 4, we show numerical curves obtained 
at different locations x in the bone. Since the BMU is created 
at time t = near x = —4.85 mm, different locations x in bone 
see the passage of the BMU at different stages of its lifetime. 
Values of x « —4.65 mm experience the passage of the BMU in 
an initiation phase with nonsteady cell densities (in particular 
a non-fully-developed osteoblast population, see footnote 11). 
Larger values of x > —3.5 mm experience the passage of 
the BMU in which cell profiles have all reached a quasi- 
steady state. 12 A possible explanation of the variability of the 
experimental data [16] reproduced in Figure 4 is therefore that 
the tetracycline data may have been collected from BMUs at 
different stages of their lifetime. This hypothesis is consistent 
with the observation of Metz et al. [16] that measured cement 
line radii exhibited great variability. If so, the frequency of 
data points lying on curves obtained from different stages of 
the lifetime of BMUs in Figure 4 may enable an experimental 
determination of the relative duration of these stages of a 
BMU's lifetime, much in the same way as the frequency of 
BMUs seen in a resorption or formation phase in histological 
cross sections enable the experimental determination of the 
duration of resorption and formation phases in BMUs [1]. 

Choosing different time intervals between the tetracycline 
injections also has a marked effect in the curves MAR versus 
R\ or R2 versus R\. This effect is shown in Figure 5. 
However, At is experimentally well controlled. By contrast, 
the knowledge of whether the BMU that formed new bone at 
the location of the cross-section at x was in an initiation stage 
or in a quasi-steady state is not straightforward to determine. 

Osteoid secretion rate. We assumed in our model that the 
osteoid secretion rate kf orm was determined by the current BMU 
cavity radius R (see Figure 2). The experimental correla- 
tion between MAR and cavity radius R does not necessarily 
imply a causal relationship between R and kf 0lm . However, 
it is possible that the local curvature of the bone substrate 
(equal to 1/R in a rotation-symmetric cortical BMU) influ- 
ences osteoblast membrane stress [20], integrin connections 
between neighbouring osteoblasts, and/or the transmission of 
osteocytic signals [21, 22]. This may influence in turn the 
osteoblasts' shape, size, and osteoid secretion rate [18,23]. 
A biological advantage of having the osteoid secretion rate 
&form effectively regulated by R is that BMU refilling becomes 
self-regulated: formation stops when refilling is complete. 
However, completion of refilling also requires that the popu- 
lation of active osteoblasts can be sustained long enough. We 
note that 'over-refilling' in cortical BMUs is not possible when 
associated with a new Haversian canal ('type I osteons') [1] 
and only possible to a certain extent when associated with 
a pre-existing Haversian canal ('type II osteons') [31,49]. 

12 An indication that the location x = —4 mm experiences the passage 
of steady-state profiles is given by the overlap of the curve obtained at 
x = —4 mm and that obtained at x = —2.5 mm in Figure 4. We note that 
the bone located at x = —2.5 mm does not experience the passage of the 
back of the BMU by t — 150 days. This explains why the curve obtained at 
x = —2.5 mm is incomplete compared to that obtained at x = —4 mm. 



Indeed, in the latter case, the Haversian canal of the old osteon 
may become smaller in the new osteon, corresponding to a net 
bone gain. 

The scarcity of experimental data on the osteoid secretion 
rate of individual osteoblasts in bone remodelling is perhaps 
surprising given that it may be easily determined by means 
of Eqs (3)-(5), combined with measurements of osteoblast 
surface densities (or secretory territory [17]) and rates of 
cavity changes. The advantage of Eq. (5) is that it holds 
locally. By measuring the secretory territory of an osteoblast 
and the (perpendicular) distance between two tetracycline 
labels at the osteoblast's position, Eq. (5) enables in principle 
the investigation of nonuniformities in the osteoid secretion 
rate kf orm going down to individual osteoblasts within the 
same BMU cross-section. The advantage of Eq. (3) is that it 
directly provides an osteoid secretion rate averaged over the 
osteoblasts in the cross-section. By measuring the surface 
area and perimeter of the BMU cavity [9], Eq. (3) enables 
the determination of osteoid secretion rates from irregular 
osteons, in constrast with Eq. (4). 

5 Conclusions 

We have developed a mathematical model of cell develop- 
ment within a cortical BMU including several influences of 
the BMU cavity on the density of active osteoblasts, either 
explicitly {e.g., via surface availability) or implicitly {e.g., via 
osteoblast-to-osteocyte transition). This model was calibrated 
against experimental data on total cell numbers and osteoblast 
surface densities at three different BMU cavity radii. The 
comparison of the calibrated model with tetracycline double 
labelling data revealed in particular the following points: 

• Tetracycline data measured on histological cross- 
sections of bone do not exhibit the start and the end 
of the refilling phase of BMUs. The build up of the 
population of active osteoblasts at the onset of refilling 
and their transition to non-synthesising bone lining cells 
at completion of refilling should be associated with low 
and nearly zero matrix apposition rates, respectively. 
The fact that this is not observed may indicate that in 
a BMU, the refilling process starts and ends rapidly, as 
suggested by our model. We note, however, that large 
and irregular Haversian canals, likely to be in an early 
stage of refilling, are often discarded from tetracycline 
analyses [12]. This difficulty could be alleviated by 
using Eq. (3), valid for any osteon shape, rather than 
Eq. (4), valid for rotation-symmetric osteons, to analyse 
the dynamics of BMU refilling. 

• The cross-sectional variability remaining in the tetra- 
cycline data reported by Metz et al. after their nor- 
malisation [16, Fig. 2] may be partly due to the fact 
that the overall refilling dynamics differs in BMUs at 
different stages of their lifetime (such as initiation stage, 
progression stage, and termination stage). Due to the 
short lifespan of osteoclasts compared to osteoblasts, the 
population of osteoclasts quickly reaches a quasi-steady 
state early in a BMU's lifetime, within a few hundreds 
of micrometres from the BMU's origination point. In 
contrast, only regions of bone about 1 mm away from 
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the BMU's origination point may experience the action of 
a fully-developed population of active osteoblasts. This 
may have two additional practical consequences, (i) 
the location of a BMU's origination may be evidenced 
by local enlargements of the Haversian pores (due to a 
local lack of osteoblasts compared to osteoclasts); (ii) it 
may be possible to determine experimentally the relative 
duration of the initiation stage and quasi-steady states of 
a BMU by counting the number of occurrences of BMUs 
whose tetracycline data falls in specific regions ("bands") 
in a plot of MAR vs cavity radius R, as suggested by 
Figure 4. 

A complete quantitative picture of the spatio-temporal dy- 
namics of BMUs based on experimental measurements remains 
to be elucidated. Whilst matrix apposition rates are often 
reported in the literature from tetracycline double labelling 
experiments, these data integrate both osteoblast number and 
osteoblast secretory activity. Experimental knowledge of 
osteoid secretion rate per osteoblast is scarce, particularly in 
relation to a BMU's internal organisation. The determination of 
cell densitites and their distribution along the different phases 
of the BMU is essential to gain insights into the biochemical 
regulation of individual cells within the BMU. The latter is 
important to fully understand how bone formation is regulated 
in bone remodelling. For example, the biochemical processes 
leading to the end of the formation phase in a BMU are not 
clear. Insights into these processes are particularly relevant 
for our understanding of osteoporosis [5], 
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Appendix A Model description 

In this appendix, we detail some aspects of the model pre- 
sented in the Methods section (Section 2). The governing 
equations of all cells and signalling molecules considered in 
the model are presented, as well as a table listing all the model 
parameters. 

Bone refilling dynamics. We first provide a derivation of 
Eq. (5), locally relating the matrix apposition rate to the 
number of osteoblasts and their level of secretory activity. An 
alternative derivation of Eq. (3) by averaging Eq. (5) over the 
BMU cross-sectional slice is also presented. 

The volume of osteoid formed at a position r of the bone 
surface during a time increment At can be characterised both 
by the geometric displacement of active osteoblasts during 
osteoid deposition, and by the amount of osteoid that was 
produced by the osteoblasts at this location during At (see 
Figure 6). The geometric displacement of an osteoblast at 
point r of the bone surface during Af is given by v 0Ba (r,t)At, 
where Vob 3 is the velocity vector of the active osteoblast, 
which describes its migration within the cross-sectional slice. 
The local volume of osteoid AVf orm (r,f) "extruded" from the 
bone surface due to this displacement is given by the volume 
of a slanted cylinder, equal to the multiplication of its base 
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Figure 6 - The displacement of the bone surface between time t and t + Af 
can be described by the displacement v B a Af of the osteoblast (OB a ) at point 
r. The volume "extruded" from the surface element dA between the times 
t and f + Af due to this displacement is equal to n ■ v Bi, AfdA. The normal 
component vi_ = n ■ Vob„ of the velocity vector corresponds to the matrix 
apposition rate. 

dA (a small, infinitesimal element of surface) and its height 
n(r,t) ■ v OBa (r, t)At, where n(r,t) is a unit vector normal to the 
bone surface (see Figure 6). On the other hand, AVf orm (r,f) is 
given by the multiplication of the volume of osteoid produced 
per osteoblast per unit time kf 0lm (r,t) and the number of active 
osteoblasts on the surface element dA. Introducing the surface 
density of active osteoblasts p OBa (r,f), one thus has: 

AVform (r,t)=n- v OB; , At dA = k foim At p 0Ba dA (15) 

The quantity = « ■ v B a describes the speed at which the 
bone surface moves in the direction normal to the surface. 
This quantity corresponds to the so-called 'matrix apposition 
rate', i.e., the thickness of the new bone layer deposited per 
unit time. From Eq. (15), the matrix apposition rate at position 
r of the bone surface and at time t is therefore equal to: 

VformC^O =feorm(»",i)PoB a (»*,?)• ( 5 ) 

To estimate the overall rate at which the BMU cavity closes 
during the refilling process, one can integrate AVf orm (7,f) in 
Eq. (15) over the bone surface of the cross-sectional slice. 
This provides the total volume of osteoid produced in the 
slice during the time increment At. This newly-formed bone 
progressively closes the BMU cavity, of surface area S(x,t). 
The change in S(x, t) is thus given by: 

AS(x,t) $ X = ~J AV form (r,f) 8xdl (16) 

where the small element of surface dA was replaced by dA = 
8x dl (dl is a length element of the cavity perimeter P). 
Assuming that kf onn , v OBj and p 0Ba do not vary significantly 
within the cross section, these quantities can be taken out 
of the integral sign when using Eq. (15), and so we retrieve 
Eq. (3) and a special case of Eqs. (5) when Af — » 0: 

§- t S(x,t) = -fc fom (x,f)p OBa (x,f)P(x,f) (3) 

= -virm(x,t)P{x,t) 

Geometric influence of surface availability. The tendency 
of cavity shrinkage to concentrate the density of osteoblasts 
is represented by the source term GOB a in the governing 
equation of active osteoblasts, Eq. (6). This effect can be 
understood intuitively as follows. During refilling, active 
osteoblasts move towards the center of the BMU as the cavity 
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radius decreases. This movement confines active osteoblasts 
into a smaller volume, and so tends to concentrate their 
density in the slice (see Figure lb, stages 2-3). An alter- 
native viewpoint is to note that the bone surface available 
to osteoblasts decreases as refilling proceeds. If the number 
of active osteoblasts on the surface is constant, this increases 
their surface density. 

The precise form of the function G(x,t) can be determined 
by considering the total number of active osteoblasts in the 
slice, 8N OBa (x,t). This total number 8N OBl , however, is 
not necessarily constant: active osteoblasts in the slice can 
be created or eliminated by biological processes. (Active 
osteoblasts are not transported into or away from the slice 
by migration from neighbouring slices since osteoblasts do 
not move along the x axis, see Section 2.1). In Eq. (6), 
biological creation or elimination of active osteoblasts in the 
slice is described in terms of changes in the volumetric density 
OB a by the source and sink terms o(x,t) = £^ b p OB p (x,f) — 

A B a OB a (x,f ) — Oqct '(x,t). At the whole slice level, the 
number of active osteoblasts created or eliminated per unit 
time is therefore described by the source and sink terms 

a(x,t)SoB a (x i t)8x since 8N 0Ba (x,t) — OB a (x,f) S OB!i (x,t) Sx, 
see Eq. (7). Hence: 



OB a 



(^ OB a)^OB a +OB a J^:SoB a 



<7 Son., Sx. 



Dividing by S OBd 8x and comparing the second equation with 
Eq. (6) gives: 



G{x,t) 



SoB a (X,t) 



R(x,t)( 



2R(x,t) 



(17) 



In Eq. (17), the second equality is obtained using Eq. (8) and 
holds only for rotation-symmetric BMUs. 

Osteoblast-to-osteocyte transition. The generation of os- 
teocytes from osteoblasts that become entrapped in the bone 
matrix during bone formation [22,52,53] induces a sink term 
~ ff ocY m tne governing equation of active osteoblasts, 
Eq. (6). Here, we determine the rate CTqcy (-M) at which 
active osteoblasts become osteocytes (in units of density 
produced per unit time) so as to reproduce a given (e.g., 
experimentally determined) osteocyte lacuna density in the 
cortical BMU, OCY exp ,(x,r), where r is the radial coordinate 
in the cross section of the (rotation-symmetric) BMU (see 
Figure lb). The total number of osteocytes in the slice at x 
at time f, /Yocy increases during bone deposition in the 
slice from zero at the onset of bone refilling t = t p (Figure lb, 
stage 1) to a constant final number at the end of bone refilling 
(Figure lb, stage 4). The correspondence between the given 
density OCY exp and the production rate of osteocytes Gq£^' 
is found by calculating Nocy(x,t) separately from OCY exp . 

and cTqcy ■ The number of osteocytes in the slice Noci(x,t) 
is given on the one hand by spatially integrating the given 
density OCY exp . from the cement line radius R c up to the cavity 
radius at time t, R(x,t) (Figure lb), and on the other hand, 
by temporally integrating the production rate <Jo°y from the 
onset of formation tp to t : 

Noc*(x,t) =2% fdrrOCY t3sp X x > r ) = /^Voct' Sqb,)(*/) 

JR(x,t) Jtp 



Differentiating the above equation with respect to t then gives: 
-2KR{x 1 t)oCY exv .{x,R(x,t))j- t R(x,t) 

= °OCY ( X J f ) ^OB a (X, t), 

whence, using Eq. (8): 

OoSfoO = f^'loB \ OCY exp.M(*>0)- (is) 

Governing equations of cell densities and signalling 
molecules concentrations. The densities of active os- 
teoblasts (OB a ), pre-osteoblast (OB p ), and active osteoclasts 
(OC a ), and the concentrations of the signalling molecules 
are governed by the material balance equation. Following 
Ref. [24], the signalling molecules considered in the model are 
transforming growth factor |3 (TGFp), parathyroid hormone 
(PTH), and the receptor-activator nuclear factor fcB system 
comprising the receptor RANK, the ligand RANKL and the de- 
coy receptor osteoprotegerin (OPG). The distribution profiles 
of uncommitted osteoblast progenitors (mesenchymal stem 
cells, OB u (x,f)) and of pre-osteoclasts (oc p (x,f)) are assumed 
to be given functions (see Table 3). Except for the dynamics 
of TGFp, which occurs on the same time scale as the resorption 
dynamics by active osteoclasts, the dynamics of the other 
signalling molecules is fast due to receptor-ligand reaction 
kinetics occurring on a much shorter time scale than the char- 
acteristic times of cell response. It can therefore be assumed 
that the concentrations of these signalling molecules quickly 
reach a pseudo-steady state (see Ref. [24] for a detailed 
discussion). This enables us to express the concentrations of 
PTH, RANKL, and OPG as algebraic equations. 

Below we summarise all the equations governing the cell 
densities and the concentrations of signalling molecules in the 
model. The material balance equation governing the dynamics 
of active osteoblasts, Eq. (6), is rewritten using the explicit 
expressions in Eqs (17) and (18), and substituting ^R(x,t) 
by the right hand side of Eq. (10). We refer the reader to 
Ref. [24] for a detailed presentation of the derivation of the 
other equations: 

J^OBafof) = % Bp (.x,f) OB p(X f ) -AoB,OB a (x,f) (19) 

+ £form(.M) OB a(*,0 ^JjOB a (jf,f) -OCY eX p. (x,R(x,t)) 

j- t OC a {x,t) = %c p (*,0 oc p(*,0-^oc a (*,0 oc a(*,f) 

-|(oc a (x,f)v OCa ), (20) 
^OBp (x, t) = @ OBu (x,t) OB u (x,t)- 3> OBp (X,t) OBp (x,t) 7 



(21) 



dt 



where 



TGFp(x,f) =n!r™ e p ^e S OC a (x,f)~£>TGFpTGFp(x,f), (22) 



-AjB„ \ X , 1 ) — i^OB„ » ( ,TGF|3 ) J 
"■OBp 

9) (rt\ — T) n-act/ RANKLfcfU 

-^OC p [X, I ) — U OC v III \ £RANKL J 7 

OCp 

^oc a (x,f)-A OCii 7r act ( i P^), 

OC a 

% Bu (x,f)=DoB u ^(^P^), 



(23) 
(24) 
(25) 
(26) 
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and 



PTH(x,f) =j3p TH /A>TH 
RANK(x,f) =A^ p NK OC p (x,f) 

/3-fOB a (,,0W^) 

£l \ K OB.rep / 

j3°|f OB a (x,f) /OFG sat +D C 



(27) 
(28) 



OPG(x,f) 



rankl(x, f ) 



x <D K 



^OB p (x,t) 



1 + kl^f- RANK(x, t ) + k^ KL OPG(x, t) 



\ K OB,act / 



(30) 



The dimensionless functions 7T act and 7T rep express the ac- 
tivation and repression of cell differentiation, apoptosis, or 
ligand expression by regulatory signalling molecules. These 
functions are determined by the fraction of receptors on the 
cell occupied by the concerned signalling molecules (see 
Refs [24,33]). Mathematically, they are "Hill" functions given 
by: 



7T re P(£) = l-7i; act £ 



1 



1+s 



(31) 



A slight change in the expression for RANKL in Eq. (30) 
has been made compared to Ref. [24]. The production of 
RANKL is now correctly proportional to the number of cells 
that express RANKL, i.e., we have replaced /Wnkl in Ref- [24, 
Eq. (19)] by j3^ NKL OB p in Eq. (30). The same inconsistency 
of having a production rate of RANKL not scaled by the number 
of osteoblasts is present in Ref. [33], which was corrected in 
Ref. [35]. The behaviour of the BMU model is not changed 
significantly by this correction. Some inconsistent behaviours 
of the temporal-only model of Ref. [33] were corrected by this 
change, see Ref. [35] for more details. 

The set of equations (19)— (3 1) together with Eq. (10) are 
solved in a co-moving frame attached to the BMU with the 
following boundary conditions: the density of pre-osteoblasts 
and active osteoblasts is zero at the (moving) front of the BMU, 
and the density of active osteoclasts is zero at the back of the 
BMU (set at a distance of 5 mm behind the front of the BMU). 
The numerical algorithm used is that of the 'method of lines' 
of Mathematica's NDSolve function. 



multiplied by a factor 500 etc. Based on Eqs (19)-(30), the 
new parameter values (denoted below by a prime) ensuring 
such scaled densities were determined by multiplying the 
previous parameter values with an adequate combination of 
the above scaling factors, according to: 





OB™ 3 *' = 


— OfoBn OB™ 3 ", 


oc max/ = 


(29) 


D 'oB u = 


«OB 

U OB u , 

«OB u 


D'oc p = 




bone' _ 
TGF p — 


"TGFp/^OCa) 


jtRANK' 

JV oc p - 


-l 


nOPG' _ 
PoB a - 


POB,, / a OBa7 


OOPG' _ 
PoBp — 




jyRANKL' 
OBp 


— iV OB p / ^OBp , 


AfRANKL' 
OB a 



^max 
T ' 



■OCa 
«OC p 

= N R 



Dr 



^OBp 3 / OfclBp, 

= AC K V«OB a , (32) 



where OB™ ax and OC™ x 



denote parameters scaling the given 
distributions of OB u s and OC p s [24]. We note that it is not 
possible to control the scaling of OB p independently from that 
of OB a without affecting their spatial relationship, and so we 
have always taken a B p = C^oBa- 



Scaling of the cell density distributions in the BMU. In 

the Results section (Section 3), the cell distribution profiles 
along x resulting from the numerical simulation are calibrated 
against experimental values of osteoblast surface densities and 
total cell numbers. To this effect, we have modified the values 
of a number of parameters compared to the parameters used 
in Ref. [24] (that were themselves taken from Ref. [33]). We 
used a scaling scheme of the model parameters such that the 
densities of OB a , OB p , OB u , OC a and OB p would scale uniformly 
without affecting much their own spatial distribution and 
without affecting their spatial relation with other cell density 
profiles. We defined scaling factors a Ba> <*oBp' <X OBu , a oc . t 
and Ctocp, such that if a B a = 500, the density of OB a s is 
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Table 3 - Model parameters and given functions 



Symbol 


Value 


Description 


OC p (x,f) 


(given function) 


pre-osteoclast density: gaussian function centred at x(t) = 
—0.35 mm + vocj, with standard deviation = 0.2 mm and 
amplitude = 2589/mm 3 . 


OB u (x,f) 


(given function) 


uncommitted osteoblast progenitors (MSC) density, equal to 

OC p (x,t) 


fcform(^) 


(given function) 


volume of osteoid formed per osteoblast per day (see Figure 2) 


^OCp 


41.26/day 


OC p — > OC a differentiation rate parameter 


^oc a 


2.82/day 


OC a apoptosis rate parameter 




81.07/day 


OB u — > OB p differentiation rate parameter 




0.166/day 


OB p — > OB a differentiation rate parameter 


A OB a 


0.0385/day 


OB a apoptosis rate 


7.RANKL 

K oc„ 

7 TGFp 

oc„ 


i nn?s . 1 n 7 mm — 3 

i.uuij iu mill 


JJalalllCLCI 1UI JvrtiNlvL, UlIlQlIlg UI1 ULp 


iiy.l mm 


parameter tor TGFp binding on OC a 


7 TGFp 
OB„ 


339 2 mm -3 

^} *j J , 111111 


narfimptpr for TfiFft riinrlincr nn OR., 

IJ CLL cllllly Lv^l 1 V / 1 1 VJ 1 jj Ulllulll 11 V / 1 1 V_/l_J|j 


7 TGFp 
"OB p 


105.6 mm" 3 


parameter for TGFp binding on OB p 


7.PTH 
"OB.act 


9.033- 10 7 mm" 3 


parameter for PTH binding on OB (for 7t' dcl ) 


J.PTH 
"OB, rep 


134039 mm" 3 


parameter for PTH binding on OB (for ;r rep ) 


7.RANKL 

"rank 


5.6655- 10" 8 mm 3 


association binding constant for RANKL and RANK 


7.RANKL 
"OPG 


1.66058- 10" 9 mm 3 


association binding constant for RANKL and OPG 


Q RANKL 
PoBp 


348. 377 /day 


production rate of RANKL per OB p 


OOPG 

Fob, 


326305/day 


production rate of OPG per OB a 


/3pTH 


1.506- 10 8 mirr 3 /day 


production rate of systemic PTH 


ivrRANKL 
JV OB p 


2.7- 10 6 


maximum number of RANKL per OB p 


ivrRANK 

iv oc p 


2326 


number of RANK receptors per OC p 


OPG sat 


1.205- 10 14 mm" 3 


OPG density at which endogeneous production stops 


£*TGFp 


0.5/day 


degradation rate of TGFp 


D RANKL 


10.13/day 


degradation rate of RANKL 


Dopg 


0.35/day 


degradation rate of OPG 


£>PTH 


86/day 


degradation rate of PTH 



«tgfp 3946 mm 3 concentration of TGFp stored in the bone matrix 

fc res 9.425 • 10~ 6 mm 3 /day volume of bone matrix resorbed per osteoclast per day 



Voc« 


30 urn/day 


longitudinal speed of active osteoclasts and of the BMU 


Rh 


20 um 


human Haversian canal radius 


Rc 


100 urn 


human cement line radius 


^OB a 


15 um 


average height of an active osteoblast 




1.25 


canine-to-human conversion factor for osteoid secretion rate 


PLC 


2300/mm 2 


surface density of bone lining cells 


OCY exp .(x,r) 


20000/mm 3 


volumetric density of osteocyte lacunae in the BMU 
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